Kevin Gao - Assignment 2: Differential Gene expression and Preliminary ORA

Overview

In the previous assignment, we took the RNASeq data from the publication Cannabidiol inhibits SARS-COV-2 replication and promotes the host innate immune response by Nguyen et al., cleaned the data, performed normalization, gene identifier mapping, and some additional preliminary analysis on the dataset. As suggested by the title, the authors of the publication studied the effects of cannabidiol on SARS-CoV-2 replication and host’s immune response. The authors hypothesized and concluded that cannabidiol inhibits SARS-CoV-2 replication by up-regulating the host IRE1α ribonuclease endoplasmic reticulum (ER) stress response and interferon signaling pathways. The experiment conditions from in the dataset can be categorized into

• Infection with SARS-CoV-2, no treatment
• Infection with SARS-CoV-2, treatment with Cannabidiol
• No infection, no treatment
• No infection, treatment with Cannabidiol

with 3 replicates for each test condition.

We obtained the dataset from GEO with ID GSE168797, associated with the study Cannabidiol inhibits SARS-COV-2 replication and promotes the host innate immune response published on Science Advances. Out of the total of 57832 genes, 13705 remained after removing low counts and genes with duplicate identifiers.

Data normalization is performed using TMM with the edgeR package. Normalization corrects the large deviation of the means of untreated groups from the groups treated with CBD, while still preserving some of the characteristics of the original sample distribution.

We observe separation of groups after normalization.

We will revisit this figure later when performing differential gene expression analysis.

Identifier mapping was performed using the package biomaRt with data from Ensembl. The Ensembl gene IDs in the original dataset are mapped to the corresponding HUGO gene symbols. The final dataset included 13705 genes with unique identifiers.

Setup

In this section, we import and install the necessary packages for this assignment, in which we will conduct a differential expression analysis using the normalized dataset and a thresholded over-representation analysis.

Differential Expression Analysis using edgeR

In this section, we are going to perform differential expression analysis.

Let us first examine the normalized data.

Normalized gene counts

	CBD_infect-1	CBD_infect-2	CBD_infect-3	Veh_infect-1	Veh_infect-2	Veh_infect-3
TSPAN6	31.887883	34.611576	36.856052	27.579725	26.797481	25.940554
DPM1	50.041084	51.301886	54.306630	45.221632	44.700001	45.116977
SCYL3	10.913215	9.376985	9.433715	18.667599	17.497179	19.244748
C1orf112	7.790084	9.467497	9.917955	17.755873	17.587255	19.335847
CFH	17.425654	19.731495	18.365254	29.699489	23.194459	29.037933
FUCA2	47.893932	54.397378	54.772936	33.460361	31.796675	31.588716
GCLC	149.679627	160.893142	162.830217	107.652101	100.186550	108.271635
NFYA	20.903687	23.514873	24.893528	36.446265	40.443930	38.421171
STPG1	16.094774	12.888830	13.253831	6.496051	6.575516	7.082978
NIPAL3	23.512212	21.306394	20.822324	12.559032	11.642267	9.770410

Sample groups

	group	lib.size	norm.factors
CBD_infect-1	CBD	673771.1	1
CBD_infect-2	CBD	677449.4	1
CBD_infect-3	CBD	685317.7	1
Veh_infect-1	Veh	782981.3	1
Veh_infect-2	Veh	781343.9	1
Veh_infect-3	Veh	793046.5	1

Workflow

We will be calculating differential expression following this general workflow:

1. Create deisng matrix. We will ensure that all factors that could contribute to differential expression is accounted for in the design matrix.

2. Use a mean-variance plot to confirm that the data follows the negative binomial distribution assumption for using the quasi-likelihood model in edgeR.

3. Fit our data to the model and calculated p-value and corrected p-value.

4. Use a threshold to extract differentially expressed genes with statistically significant differences.

5. Plot and visualize.

Creating the Design Matrix

In order to perform statistical testing, we need a design matrix the defines our model. Notice that in our dataset, there are three factors:

1. Treatment with CBD
2. Infection status (infected or not)
3. Patient

Hence, ideally, we would like to account for all three factors in our design matrix.

Design matrix

(Intercept)	samples$patient2	samples$patient3	samples$cbd_treatmentCBD	samples$infectedmock
1	0	0	1	0
1	1	0	1	0
1	0	1	1	0
1	0	0	0	0
1	1	0	0	0
1	0	1	0	0
1	0	0	1	1
1	1	0	1	1
1	0	1	1	1
1	0	0	0	1
1	1	0	0	1
1	0	1	0	1

Distribution of Data

For our downstream analysis, we are going to use edgeR. We chose edgeR because it is specifically designed for RNASeq data. However, one important underlying assumption for using the quasi-likelihood model is that the data follows a negative binomial distribution. We need to verify that our dataset indeed meets that assumption.

Let us the calculate the dispersion and plot to visualize the mean-variance relationship.

edgeR::plotMeanVar(d,
                   show.raw.vars = TRUE,
                   show.tagwise.vars = TRUE,
                   NBline = TRUE,
                   show.ave.raw.vars = TRUE,
                   show.binned.common.disp.vars = TRUE,
                   main = "Mean-Variance Plot for Variance and Dispersion of Our Data")
# display legend
legend("topleft", 
       legend=c("Raw Data", "Tagwise Dispersion", "Average Raw Variances", 
                "Binned Common Dispersion", "Negative Binomial Line"), 
       col = c("grey", "lightblue", "maroon", "red", "dodgerblue2"), pch=c(1,1,4,4,NA), lty=c(0,0,0,0,1), lwd=c(1,1,1,1,2), cex=0.6)

As demonstrated by the mean-variance plot above, we can see that the dispersion and variance of our data indeed roughly follows the negative binomial distribution.

Analysis Using edgeR

Now, we have created the design matrix and verified the assumption for the data to be negative-binomially distributed, we can proceed to the next stage of our analysis and perform statistical testing and corrections to ensure that we only get significantly differentially expressed genes. We used the quasi-likelihood models since our dataset is from an RNASeq experiment and quasi-likelihood models are best suited to handle RNASeq data.

Once we have fit the model, we can proceed to calculate differential expression. Recall that our goal is to verify the role of cannabidiol in affecting the replication ability of SARS-CoV-2, so we will be using cbd_treatment as the contrast.

Top differentially expressed genes

	logFC	logCPM	F	PValue	FDR
DNAJB9	3.296477	6.632298	2376.139	1.174076e-13	1.225637e-09
IGFBP1	4.315555	7.225448	2191.283	1.788598e-13	1.225637e-09
HSPA5	3.721184	11.833327	1920.096	3.552926e-13	1.623095e-09
CRELD2	3.023632	6.669438	1395.860	1.858463e-12	5.742205e-09
KIF20A	-2.327554	5.891384	1344.930	2.253309e-12	5.742205e-09
ADGRF4	2.315894	4.618755	1287.278	2.827571e-12	5.742205e-09
GABARAPL1	2.351451	6.414219	1278.222	2.932903e-12	5.742205e-09
NDRG4	1.776359	6.129025	1133.516	5.463712e-12	9.159587e-09
HERPUD1	3.400467	8.195296	1061.159	7.686235e-12	9.159587e-09
STC2	2.742212	6.444858	1059.247	7.758262e-12	9.159587e-09

We can examine the number of genes pass the threshold and correction. We are using 0.05 as the threshold for p-value as it is commonly used in practice.

The threshold of 0.05 gives us quite a lot of genes, in order to get more meaningful hits in the downstream gene enrichment analysis, we will make the thershold more stringent.

Visualization of Differentially Expressed Genes

Now we can retrieve the list of differential expressed genes and visualize using different plots. We will first plot them on a volcano plot. Each gene is represented by a point in the plot. The horizontal axis of the plot is the \(\log_2\) fold change and the vertical axis is the \(-\log_{10}p\) which indicates the statistical significance of each gene (how likely the differential is due to actual biological variation).

Next, we will visualize the upregulated and downregulated genes across different test conditions using a heatmap.

Discussion

1. Initially, I used the p-value of 0.05 as it is widely used in practice. This gives us 5501 genes prior to correction. This is quite a large number of genes, so we changed the p-value threshold to 0.01 to limit the number of genes included. We further added the criteria that a gene must have an absolute log fold change greater than 1.5. By doing so, we are contraining ourselves to only get the genes that are highly differentially expressed with high statistical significance.

2. For correction, we used Benjamini-Hochberg. The two main methods discussed for correcting family-wise false discovery rate are Bonferroni and Benjamini-Hochberg correction. We want to get meaningful hits without excluding statistically significant ones. Bonferroni’s method is overly stringent so it is not quite suitable for our purposes. Using Benjamini-Hochberg will give us a richer set of genes that we can use for downstream analysis.

3. Volcano plot shown above. Note that we also applied the restriction that only genes with aboslute log fold change greater than 1.5 are to be included.

4. There is significant clustering within conditions. This suggests that the test condition (CBD v.s. no CBD) does have a significant effect on gene expressions.

Thresholded Overrepresentation Analysis using g:Profiler

For the final part of this assignment, we will perform a thresholded overrepresentation analysis using g:Profiler. In the previous section, we have compiled a list of differentially expressed genes. Here, we want to further divide them into upregulated and downregulated genes.

upregulated_gene_lst <- diff_exp_lst[diff_exp_lst$logFC > 0,]
downregulated_gene_lst <- diff_exp_lst[diff_exp_lst$logFC < 0,]

We use the R package for g:Profiler to perform the gene enrichment analysis. For correction, we used FDR as it is less stringent than Bonferroni and is introduced as the preferred correction method in class. We used GO Biological Process, GO Molecular Function, and KEGG as those are the ones used by the author of the original publication. For a more detailed overview of the workflow, please refer to the Discussion subsection located at the end of this section.

Upregulated Genes

Top genesets using list of upregulated genes

term_name	term_id	source
response to endoplasmic reticulum stress	GO:0034976	GO:BP
proteasomal protein catabolic process	GO:0010498	GO:BP
ERAD pathway	GO:0036503	GO:BP
ubiquitin-dependent ERAD pathway	GO:0030433	GO:BP
proteasome-mediated ubiquitin-dependent protein catabolic process	GO:0043161	GO:BP
autophagy	GO:0006914	GO:BP
process utilizing autophagic mechanism	GO:0061919	GO:BP
response to topologically incorrect protein	GO:0035966	GO:BP
macroautophagy	GO:0016236	GO:BP
endoplasmic reticulum to Golgi vesicle-mediated transport	GO:0006888	GO:BP

For context, let’s examine the top term from each data source.

Top terms from each data source using list of upregulated genes

	term_name	term_id	source
1	response to endoplasmic reticulum stress	GO:0034976	GO:BP
331	unfolded protein binding	GO:0051082	GO:MF
401	Protein processing in endoplasmic reticulum	KEGG:04141	KEGG

We can visualize the distribution of top terms from each data source using an Manhattan plot. This is the distribution of terms prior to removing large terms with over 500 genes.

gprofiler2::gostplot(up_top_terms_all) %>% plotly::layout(title = "Manhattan plot showing distribution of terms \nfrom each data source using list of upregulated genes", font = list(size = 10))

Downregulated Genes

We do the same for the downregualted genes.

Top genesets using list of downregulated genes

term_name	term_id	source
mitotic nuclear division	GO:0140014	GO:BP
sister chromatid segregation	GO:0000819	GO:BP
mitotic sister chromatid segregation	GO:0000070	GO:BP
nuclear division	GO:0000280	GO:BP
chromosome segregation	GO:0007059	GO:BP
organelle fission	GO:0048285	GO:BP
nuclear chromosome segregation	GO:0098813	GO:BP
cell cycle phase transition	GO:0044770	GO:BP
regulation of cell cycle phase transition	GO:1901987	GO:BP
DNA replication	GO:0006260	GO:BP

Top terms from each data source using list of downregulated genes

	term_name	term_id	source
1	mitotic nuclear division	GO:0140014	GO:BP
474	chromatin binding	GO:0003682	GO:MF
535	Cell cycle	KEGG:04110	KEGG

Plot the Manhattan plot showing distribution of terms from each data source using list of downregulated genes.

gprofiler2::gostplot(down_top_terms_all) %>% plotly::layout(title = "Manhattan plot showing distribution of terms\n from each data source using list of downregulated genes", font = list(size = 10))

All Differentially Expressed Genes

Finally, for all differentially expressed genes.

Top genesets using list of all differentially expressed genes

term_name	term_id	source
mitotic sister chromatid segregation	GO:0000070	GO:BP
mitotic nuclear division	GO:0140014	GO:BP
sister chromatid segregation	GO:0000819	GO:BP
response to endoplasmic reticulum stress	GO:0034976	GO:BP
proteasomal protein catabolic process	GO:0010498	GO:BP
organelle fission	GO:0048285	GO:BP
chromosome segregation	GO:0007059	GO:BP
nuclear chromosome segregation	GO:0098813	GO:BP
nuclear division	GO:0000280	GO:BP
regulation of cell cycle phase transition	GO:1901987	GO:BP

Top terms from each data source using list of all differentially expressed genes

	term_name	term_id	source
1	mitotic sister chromatid segregation	GO:0000070	GO:BP
514	cadherin binding	GO:0045296	GO:MF
583	Protein processing in endoplasmic reticulum	KEGG:04141	KEGG

gprofiler2::gostplot(top_terms_all) %>% plotly::layout(title = "Manhattan plot showing distribution of terms from each data source", font = list(size = 10))

Discussion

1. We used g:Profiler as we have extensively discussed about it in class. It has a nice web-based interface as well as easy-to-use APIs in the form of an R package. Furthermore, the data sources on g:Profiler is frequently updated and include the data sources that we are most interested in.

2. We used GO Biological Process, GO Molecular Function, and KEGO. We chose these data sources since they are also used by the author of the original paper from which I obtained the dataset. However, since KEGO is a commercial data source, if it was not for that fact that the original paper used it, I would personally refrain from using it. The original author also used other data sources including Canonical Pathways but since it is not part of g:Profiler, we did not include these data sources. The annotation source versions are as follows: Ensembl 105, Ensembl Genomes 52 (database built on 2022-02-14)

3. For all three analysis (using upregulated, downregulated, all differentially expressed), we used a threshold between 2 and 500. We set the upper bound to 500 because we do not want to include overly broad and generic terms that will not give us meaningful insights into the roles of the differentially expressed genes. The set of upregulated genes returned 433 genesets; the set of downregulated genes returned 560 genesets; the set of all differentially expressed genes returned 616 gene sets.

4. The result using the whole list (set of all differentially expressed genes) is more predominantly represented by the downregulated genes. It also does not provide a lot of insights into the roles of the genes as it is a mix of drastically different and seemingly unrelated terms (e.g. mitotic nuclear division, organelle fission, response to endoplasmic reticulum stress, etc.). However, by using separating into upregualted and downregulated genes and performing the gene enrichment analysis separately, we get more meaningful terms: the upregualted genes are mostly associated with pathways and processes involving endoplasmic reticulum whereas the downregulated genes are mostly associated with the cell cycle.

Interpretation

1. Yes. The overrepresentation analysis results support the conclusion discussed in the original paper. The original paper claims that CBD treatment can inhibit SARS-CoV-2 replication through the host’s ER stress response pathway. Although we were not able to verify whether or not this has any affect on SARS-CoV-2 replication, we are able to confirm based on the top terms returned by the overrepresentation analysis that the genes associated with ER stress response are indeed upregulated in the samples treated with CBD.

2. The conclusion that CBD is correlated to reduced risk or severity of SARS-CoV-2 is supported by testing adult patients as presented by the original paper. There are also other papers who discussed the role of CBD in SARS-CoV-2 infection (van Breeman et al.). The result from our overrepresentation analysis is supported by the original paper, in which the author has also performed gene enrichment analysis with similar outcomes. However, since this is still an early research in this topic and SARS-CoV-2 continues to evolve, there does not seem to be other unrelated groups who could provide additional supporting evidence to our results from the overrepresentation analysis besides the original paper.

Journal

Link to my journal entry for this assignment: https://github.com/bcb420-2022/Kevin_Gao/wiki/Assignment-2

References

• Chen Y, Lun ATL, Smyth GK (2016). From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research 5, 1438
• Durinck, S., Spellman, P. T., Birney, E., & Huber, W. (2009). Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nature - protocols, 4(8), 1184–1191. https://doi.org/10.1038/nprot.2009.97
• Martin Morgan (2021). BiocManager: Access the Bioconductor Project Package Repository. R package version 1.30.16. https://CRAN.R-project.org/package=BiocManager
• McCarthy DJ, Chen Y and Smyth GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research 40, 4288-4297
• Nguyen, L. C., Yang, D., Nicolaescu, V., Best, T. J., Ohtsuki, T., Chen, S.-N., Friesen, J. B., Drayman, N., Mohamed, A., Dann, C., Silva, D., Gula, H., Jones, K. A., Millis, J. M., Dickinson, B. C., Tay, S., Oakes, S. A., Pauli, G. F., Meltzer, D. O., … Rosner, M. R. (2021). Cannabidiol inhibits SARS-COV-2 replication and promotes the host innate immune response. Science Advances. https://www.science.org/doi/abs/10.1126/sciadv.abi6110
• Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-140
• Gu, Z. (2016) Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics.
• Kolberg L, Raudvere U, Kuzmin I, Vilo J, Peterson H (2020). “gprofiler2- an R package for gene list functional enrichment analysis and namespace conversion toolset g:Profiler.” F1000Research, 9 (ELIXIR)(709). R package version 0.2.1.
• Stefan Milton Bache and Hadley Wickham (2020). magrittr: A Forward-Pipe Operator for R. https://magrittr.tidyverse.org, https://github.com/tidyverse/magrittr.
• van Breemen, R. B., Muchiri, R. N., Bates, T. A., Weinstein, J. B., Leier, H. C., Farley, S., & Tafesse, F. G. (2022). Cannabinoids Block Cellular Entry of SARS-CoV-2 and the Emerging Variants. Journal of natural products, 85(1), 176–184. https://doi.org/10.1021/acs.jnatprod.1c00946